logo of company

Intersections of Adversity & Neurodiversity: Adverse Childhood Experiences’ Association to Mental Health & the Buffering Role of Flourishing

This study extends analyses conducted by Bethell et al. (2023), who explored the role of family resilience and connection in promoting flourishing among U.S. children facing adversity. Our study applies a similar framework but focuses on neurodivergent children, examining how adverse childhood experiences (ACEs) influence mental health outcomes and how flourishing behaviors—such as curiosity, emotional control, and task persistence—may buffer these effects. This study was conducted by the Slopen Laboratory at the Harvard T.H. Chan School of Public Health.

Author: Adrián Medina

Date: October 7, 2024

Project Overview

This repository contains files relevant to a study investigating the relationship between adverse childhood experiences (ACEs) and mental health outcomes in children with neurodevelopmental disorders (NDDs), such as Autism Spectrum Disorder (ASD) and ADHD. The study specifically examines the moderating role of child flourishing behaviors—curiosity in learning, emotional control, and task persistence—on the effects of ACEs on mental health outcomes.

Study Design:
The study utilizes data from the National Survey of Children’s Health (NSCH), focusing on children aged 6-17 with reported neurodevelopmental issues (N = 44,776, MAge = 12.2). This cross-sectional analysis examines ACE exposure and its impact on three primary mental health outcomes: anxiety, depression, and behavioral issues. The moderating role of child flourishing behaviors is assessed using interaction terms in logistic regression models.

Data Collection:
Key variables collected include:

  • Adverse Childhood Experiences (ACEs): Derived from parent-reported experiences such as exposure to violence, household dysfunction, and discrimination.
  • Mental Health Outcomes: Parent-reported diagnoses of anxiety, depression, and behavioral problems, confirmed by health care providers or educators.
  • Child Flourishing Metrics: Based on parent-reported behaviors such as curiosity in learning, emotional control, and task persistence, which were used to develop a Child Flourishing Index (CFI).

Analytic Approach:
The primary analytic approach involved logistic regression models to assess the association between ACEs and mental health outcomes, with interaction terms to explore the moderating effect of child flourishing. Key covariates include age, sex, race/ethnicity, and socioeconomic status.

Goals:
The study aims to:

  • Examine the dose-response relationship between the number of ACEs and risks of mental health challenges in neurodivergent children.
  • Identify the protective role of child flourishing behaviors, offering insights into resilience-promoting interventions for this vulnerable population.
Tip

For detailed information on the dataset variables, refer to the NSCH Data Dictionary.

Code Workflow

  1. Data Frame Initialization
  2. Analytic Data Preparation
  3. Frequencies & Descriptive Statistics
  4. Inferential Statistics
  5. Session Information

Data Frame Initialization

Set up the R environment by configuring the CRAN repository, installing essential packages, setting the base path, and loading the primary dataset into the local environment.

Expand Code 
# Set CRAN repository for consistent package installation
options(repos = c(CRAN = "http://cran.rstudio.com/"))

# Install and load the pacman package for efficient package management
if (!require(pacman)) install.packages("pacman")
library(pacman)

# Use p_load function to install (if necessary) and load packages
p_load(dplyr, tidyr, tidyverse, knitr, ggplot2, Hmisc, broom, stats, kableExtra, webshot2, sjPlot, margins, gtsummary, ggstats, broom.helpers, patchwork, sjlabelled, sjmisc)

# Specify the 'base_path' where you can find your data files, ASSUMING they're in the same directory, & set it as WD
base_path <- "~/GitHub/Adversity-Neurodiversity/data_files"
setwd(base_path)

# Load primary data file
neurodivergent_data <- read_csv("neurodivergent_data.csv")
1
See details under the Data Extraction & Filtering (Archived Code) callout below!
Data Extraction & Filtering (Archived Code)

The file being used for these analyses is a subset of a “master” file, Data2_2016to2022.csv, containing the compiled contents of the data releases from the National Survey on Child Health’s 2016-2022 cycles. The master file is nearly 500 MB (0.5 GB) in size, making it cumbersome to process in real-time. For efficiency and ease of analysis, I am using a smaller subset that contains only the relevant data. This subset was generated using the data extraction and filtering process detailed in the archived code below, in an effort to maintain transparency and ensure reproducibility of the workflow.

Expand Code 
# Define the path to the dataset
data_path <- "~/Downloads/Data2_2016to2022.csv"  
Data2_2016to2022 <- read.csv(data_path)

# Select specific variables, create 'neurodivergent' variable, & filter data based on age and neurodivergent status
neurodivergent_data <- Data2_2016to2022 %>%
  select(year, fpl, fpl_mean, sc_age_years, sc_hispanic_r, sc_race_r, sc_sex, higrade_tvis, k2q35a, k2q35c, k2q38a,
         k2q38c, k2q30a, k2q30c, k2q36a, k2q36c, k2q60a, k2q60c, k2q37a, k2q37c, downsyn, downsyn_desc, 
         k2q31a, k2q31c, k2q61a, cerpals_desc, k2q33a, k2q33b, k2q33c, k2q32a, k2q32b, k2q32c, k2q34a, k2q34b, k2q34c, 
         ace1, ace3, ace4, ace5, ace6, ace7, ace8, ace9, ace10, ace11, ace12, k6q71_r, k7q85_r, k7q84_r) %>%
  # Creating 'neurodivergent' variable as NSCH does not explicitly inquire about neurodivergent status.
  # Using a list of diagnoses provided by NSCH to define neurodivergent individuals.
  mutate(neurodivergent = if_else(k2q35a == 1 | k2q38a == 1 | k2q36a == 1 | k2q60a == 1 | 
                                  k2q37a == 1 | k2q30a == 1 | downsyn == 1 | 
                                  k2q31a == 1 | k2q61a == 1, 1, 0)) %>%
  # Filter data to include only individuals aged 6 or above and identified as neurodivergent
  filter(sc_age_years >= 6, neurodivergent == 1)

# Resulting filtered data to be used for further analysis
write.csv(neurodivergent_data, "neurodivergent_data.csv")

Analytic Data Preparation

Set up recoding of variables by transforming ‘ace1’ into a dichotomous variable, calculating the Childhood Flourishing Index (CFI), and categorizing total ACE counts. Recode education, income, sex, race/ethnicity, and CFI into categorical factors for analysis.

Expand Code 
# Recode 'ace1' to dichotomous variable based on analytic specifications
neurodivergent_data$ace1_recode <- ifelse(neurodivergent_data$ace1 %in% c(2, 3, 4), 1, 
                                          ifelse(neurodivergent_data$ace1 == 1, 0, NA))

# Recode ACE counts and calculate Childhood Flourishing Index (CFI)
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    # Total ACE count for each child excluding missing values
    ace_total = rowSums(select(., ace1_recode, ace3:ace12) == 1, na.rm = TRUE),
    
    # Categorize total ACEs for further analysis
    ace_counts = cut(ace_total, breaks = c(-1, 0, 1, 3, Inf), labels = c('0', '1', '2-3', '4+'), right = TRUE),
    ace_counts = factor(ace_counts, levels = c("0", "1", "2-3", "4+")), # Recreate 'ace_counts' factor

    # Calculate the Childhood Flourishing Index (CFI)
    CFI = rowSums(select(., k6q71_r, k7q85_r, k7q84_r) == 1, na.rm = TRUE),
    
    # Recode CFI into dichotomous variable
    CFI_dich = case_when(
      CFI == 3 ~ "Flourishing",
      CFI %in% 0:2 ~ "Not Flourishing",
      TRUE ~ NA_character_
    ),
    CFI_dich = factor(CFI_dich, levels = c("Not Flourishing", "Flourishing")) # Recreate 'CFI_dich' factor
  )

# Recode binary outcomes for anxiety, depression, and behavioral issues into factors
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    # Recode Anxiety (k2q33b)
    Anxiety = factor(k2q33b, levels = c(2, 1), labels = c("No", "Yes")),
    
    # Recode Depression (k2q32b)
    Depression = factor(k2q32b, levels = c(2, 1), labels = c("No", "Yes")),
    
    # Recode Behavioral Issues (k2q34b)
    Behavioral_Issues = factor(k2q34b, levels = c(2, 1), labels = c("No", "Yes"))
  )

# Recode education and income categories
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    # Recode education into categorical factors
    highgrade_tvis_cat = case_when(
      higrade_tvis == 1 ~ "Less than high school",
      higrade_tvis == 2 ~ "High school (including vocational)",
      higrade_tvis == 3 ~ "Some college or associate degree",
      higrade_tvis == 4 ~ "College degree or higher",
      TRUE ~ NA_character_
    ),
    highgrade_tvis_cat = factor(highgrade_tvis_cat, levels = c("Less than high school", "High school (including vocational)", "Some college or associate degree", "College degree or higher")),

    # Recode federal poverty level categories into factors
    fpl_cat = case_when(
      fpl %in% c(50:99) | fpl_mean < 100 ~ "Less than 100%",
      fpl %in% c(100:199) | fpl_mean < 200 ~ "100% to 199%",
      fpl %in% c(200:299) | fpl_mean < 300 ~ "200% to 299%",
      fpl %in% c(300:399) | fpl_mean < 400 ~ "300% to 399%",
      fpl %in% c(400:999) | fpl_mean < 999 ~ "400% or greater",
      TRUE ~ NA_character_
    ),
    fpl_cat = factor(fpl_cat, levels = c("Less than 100%", "100% to 199%", "200% to 299%", "300% to 399%", "400% or greater"))
  )

# Recode sex and race categories
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    # Recode sex into categorical factors
    sc_sex_cat = case_when(
      sc_sex == 1 ~ "Male",
      sc_sex == 2 ~ "Female",
      TRUE ~ NA_character_
    ),
    sc_sex_cat = factor(sc_sex_cat, levels = c("Male", "Female")),

    # Recode race and ethnicity
    race = case_when(
      sc_race_r == 1 ~ "White",
      sc_race_r == 2 ~ "Black",
      sc_race_r == 3 ~ "American Indian or Alaska Native",
      sc_race_r %in% 4:5 ~ "Asian, Native Hawaiian, or Pacific Islander",
      sc_race_r %in% 6:7 ~ "Other or mixed race",
      TRUE ~ NA_character_
    ),

    ethnicity = case_when(
      sc_hispanic_r == 1 ~ "Hispanic/Latino",
      sc_hispanic_r == 2 ~ "Non-Hispanic/Latino",
      TRUE ~ NA_character_
    ),

    # Recreate 'sc_race_cat' from race and ethnicity combinations
    sc_race_cat = case_when(
      race %in% c("White") & ethnicity == "Non-Hispanic/Latino" ~ "Non-Hispanic White",
      race %in% c("Black") & ethnicity == "Non-Hispanic/Latino" ~ "Non-Hispanic Black or African American",
      race %in% c("American Indian or Alaska Native") & ethnicity == "Non-Hispanic/Latino" ~ "Non-Hispanic American Indian or Alaska Native",
      race %in% c("Asian, Native Hawaiian, or Pacific Islander") & ethnicity == "Non-Hispanic/Latino" ~ "Non-Hispanic Asian, Native Hawaiian, or Pacific Islander",
      race %in% c("Other or mixed race") & ethnicity == "Non-Hispanic/Latino" ~ "Non-Hispanic Other or mixed race",
      ethnicity == "Hispanic/Latino" ~ "Hispanic or Latino of any race",
      TRUE ~ NA_character_
    ),
    sc_race_cat = factor(sc_race_cat, levels = c("Non-Hispanic White", "Non-Hispanic Black or African American", 
                                                "Non-Hispanic American Indian or Alaska Native", "Non-Hispanic Asian, Native Hawaiian, or Pacific Islander", "Non-Hispanic Other or mixed race", "Hispanic or Latino of any race"))
  )
1
ace1 is originally coded by NSCH as a frequency measure, but it’s needed as indicative variable of presence like the other ace# variables.
2
Part of this includes the calculating of the Childhood Flourishing Index (CFI) based on Bethell et al. (2019) criteria

Frequencies & Descriptive Statistics

Expand Code 
# Constructing the frequency table for Anxiety, Depression, and Behavioral Issues where the response is "Yes"
neurodiv_freq_table <- neurodivergent_data %>%
  dplyr::summarise(
    Anxiety_Yes = sum(Anxiety == "Yes", na.rm = TRUE),
    Depression_Yes = sum(Depression == "Yes", na.rm = TRUE),
    Behavioral_Issues_Yes = sum(Behavioral_Issues == "Yes", na.rm = TRUE)
  )

# Display table in a simple HTML format using kable
neurodiv_freq_table %>%
  knitr::kable("html", col.names = c("Anxiety", "Depression", "Behavioral Issues"), 
               caption = "Frequency of Mental Health Conditions in Neurodivergent Sample")
Frequency of Mental Health Conditions in Neurodivergent Sample
Anxiety Depression Behavioral Issues
13725 6289 13535
Expand Code 
# Count the number of subjects with a CFI score of 3 using a logical condition, as this is considered 'Flourishing.'
num_subjects_cfi_3 <- sum(neurodivergent_data$CFI == 3, na.rm = TRUE)
cat("Number of Subjects Flourishing (CFI Score of 3):", num_subjects_cfi_3, "\n")
Number of Subjects Flourishing (CFI Score of 3): 3265
Expand Code 
# Create the summary avg risk 'Y' data frame
neurodiv_age_summary <- neurodivergent_data %>%
  dplyr::summarise(
    score_mean = mean(sc_age_years, na.rm = TRUE),
    n = n(), # Sample size for each group
    sem = sd(sc_age_years, na.rm = TRUE) / sqrt(n()),
    .groups = 'drop'
  )

# Raincloud plot of ages for neurodivergent participants with color
ggplot(neurodivergent_data, aes(x = 1, y = sc_age_years)) + # Set x to 1, as there's only one group
  PupillometryR::geom_flat_violin(adjust = 1.5, trim = FALSE, alpha = 0.5, color = NA, fill = "deeppink", position = position_nudge(x = 0.1, y = 0)) +
  geom_point(alpha = 0.1, position = position_jitter(width = .05), size = .10, shape = 20, color = "deeppink") +
  geom_boxplot(outlier.shape = NA, alpha = 0, width = 0.1, position = position_dodge(width = 0.3), colour = "black") +
  geom_point(data = neurodiv_age_summary, aes(x = 1, y = score_mean), shape = 18, color = "deeppink", size = 3, position = position_nudge(x = 0.1)) +
  geom_errorbar(data = neurodiv_age_summary, aes(x = 1, y = score_mean, ymin = score_mean - sem, ymax = score_mean + sem), width = 0.05, color = "deeppink", position = position_nudge(x = 0.1)) +
  labs(
    title = "Age Distribution of Neurodivergent Participants", 
    y = "Age (Years)", 
    x = "") + # No x-axis label since there's only one group
  scale_y_continuous(breaks = seq(5, 18, by = 1), limits = c(5, 18)) + # Set y-axis range and increments
  theme_minimal() +
  theme(axis.text.x = element_blank()) # Remove x-axis text

Expand Code 
# Bar plot of ACE counts categories
ggplot(neurodivergent_data, aes(x = factor(ace_counts), fill = factor(ace_counts))) +
  geom_bar(show.legend = FALSE) +
  geom_text(stat = 'count', aes(label = ..count..), position = position_stack(vjust = 0.5), color = "white") +
  labs(x = "ACE Counts", y = "Number of Children", fill = "ACE Counts",
       title = "Distribution of Children by ACE Counts") +
  scale_x_discrete(labels = c('0' = "0 ACEs", '1' = "1 ACE", '2-3' = "2-3 ACEs", '4+' = "4+ ACEs")) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 0))

Expand Code 
# Bar plot of CFI counts categories
ggplot(neurodivergent_data, aes(x = factor(CFI), fill = factor(CFI))) +
  geom_bar(show.legend = FALSE) +
  geom_text(stat = 'count', aes(label = ..count..), position = position_stack(vjust = 0.5), color = "white") +
  labs(x = "Childhood Flourishing Index (CFI) Score",
       y = "Number of Children",
       fill = "CFI Score",
       title = "Distribution of Children by CFI Score") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 0))

Expand Code 
# Bar plot of CFI_dich counts categories
ggplot(neurodivergent_data, aes(x = CFI_dich)) +
  geom_bar(aes(fill = CFI_dich), show.legend = FALSE) +
  geom_text(stat = 'count', aes(label = ..count..), position = position_stack(vjust = 0.5), color = "white") +
  labs(title = "Distribution of CFI Dichotomous Variable",
       x = "Child Flourishing Index (Dichotomous)",
       y = "Count",
       fill = "CFI") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 0))

Expand Code 
# Bar plot for educational attainment counts categories
ggplot(neurodivergent_data %>% filter(!is.na(highgrade_tvis_cat)), aes(x = highgrade_tvis_cat)) +
  geom_bar(aes(fill = highgrade_tvis_cat), show.legend = FALSE) +
  geom_text(stat = 'count', aes(label = ..count..), position = position_stack(vjust = 0.5), color = "white") +
  labs(title = "Distribution of Adults' Highest Educational Attainment",
       x = "Educational Attainment",
       y = "Count",
       fill = "Educational Attainment") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Expand Code 
# Bar plot for sex counts distribution
ggplot(neurodivergent_data, aes(x = sc_sex_cat)) +
  geom_bar(aes(fill = sc_sex_cat), show.legend = FALSE) +
  geom_text(stat = 'count', aes(label = ..count..), position = position_stack(vjust = 0.5), color = "white") +
  scale_fill_manual(values = c("Male" = "lightblue", "Female" = "pink")) +
  labs(title = "Distribution of Sex",
       x = "Sex",
       y = "Count") +
  theme_minimal()

Expand Code 
#Bar plot for SES counts categories
ggplot(neurodivergent_data, aes(x = fpl_cat)) +
  geom_bar(aes(fill = fpl_cat), show.legend = FALSE) +
  geom_text(stat = 'count', aes(label = ..count..), position = position_stack(vjust = 0.5), color = "white") +
  labs(title = "Distribution of Socioeconomic Status",
       x = "Socioeconomic Status",
       y = "Count") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Expand Code 
# Bar plot for race counts categories
ggplot(neurodivergent_data, aes(x = sc_race_cat)) +
  geom_bar(aes(fill = sc_race_cat), show.legend = FALSE) +
  geom_text(stat = 'count', aes(label = ..count..), vjust=-0.3, size=3.5) +
  labs(title = "Distribution by Race",
       x = "Race",
       y = "Count") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Expand Code 
# Calculate mean and standard deviation of age for neurodivergent participants
age_stats <- neurodivergent_data %>%
  filter(neurodivergent == 1) %>%
  summarise(
    Mean_Age = mean(sc_age_years, na.rm = TRUE),
    SD_Age = sd(sc_age_years, na.rm = TRUE)
  )

# Print the statistics using a nicely formatted table
knitr::kable(age_stats, caption = "Age Statistics of Neurodivergent Participants")

# Calculate prevalence (%) of diagnoses among subject sample
variables <- c("k2q31a", "k2q35a", "k2q38a", "k2q30a", "k2q36a", "k2q60a", "k2q37a", "downsyn", "k2q61a")
variable_names <- c("ADHD", "Autism/ASD", "Tourette’s Syndrome", "Learning Disability", 
                    "Development Delay", "Intellectual Disability", 
                    "Speech/Other Language Disorder", "Down Syndrome", "Cerebral Palsy")

# Pre-calculate total respondents to avoid repeated computation
total_respondents <- colSums(!is.na(neurodivergent_data[variables]))

# Vectorized prevalence calculation
prevalences <- colSums(neurodivergent_data[variables] == 1, na.rm = TRUE) / total_respondents * 100
prevalence_table <- data.frame(Diagnosis = variable_names, Prevalence = prevalences)

# Display prevalence data in a table
knitr::kable(prevalence_table, caption = "Prevalence of Neurodivergent Conditions")

# Prepare the data for cross-tabulations
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    k2q32b = factor(k2q32b, levels = c(2, 1), labels = c("No", "Yes")),
    k2q33b = factor(k2q33b, levels = c(2, 1), labels = c("No", "Yes")),
    k2q34b = factor(k2q34b, levels = c(2, 1), labels = c("No", "Yes")),
    CFI_dich = factor(CFI_dich, levels = c(0, 1), labels = c("Not flourishing", "Flourishing")),
    ace_counts = factor(ace_counts, levels = c(1, 2, 3, 4), labels = c("0", "1", "2 or 3", "4+"))
  )

# Create the cross-tabulations
table_data <- neurodivergent_data %>%
  select(k2q32b, k2q33b, k2q34b, CFI_dich, ace_counts) %>%
  pivot_longer(cols = c(k2q32b, k2q33b, k2q34b, CFI_dich), names_to = "Condition", values_to = "Status") %>%
  group_by(Condition, Status, ace_counts) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(Condition, ace_counts) %>%
  mutate(percent = n / sum(n) * 100) %>%
  pivot_wider(names_from = ace_counts, values_from = c(n, percent), names_sep = ", ")

# Format the table for presentation using 'kable' and 'kableExtra'
knitr::kable(table_data, format = "latex", caption = "Frequency and Percentage of Psychopathology and Flourishing across ACEs Categories") %>%
  kableExtra::kable_styling("striped", full_width = FALSE)

Inferential Statistics

Expand Code 
# Logistic regression models without interactions
model_anx <- glm(Anxiety ~ ace_counts + sc_race_cat + sc_sex_cat + fpl_cat + highgrade_tvis_cat + sc_age_years,
               data = neurodivergent_data, family = binomial)

model_dep <- glm(Depression ~ ace_counts + sc_race_cat + sc_sex_cat + fpl_cat + highgrade_tvis_cat + sc_age_years,
               data = neurodivergent_data, family = binomial)

model_beh <- glm(Behavioral_Issues ~ ace_counts + sc_race_cat + sc_sex_cat + fpl_cat + highgrade_tvis_cat + sc_age_years,
               data = neurodivergent_data, family = binomial)

# Create a summary table for the Anixety model
tbl_anx <- tbl_regression(model_anx, exponentiate = TRUE, label = list(
      ace_counts ~ "ACEs",
      fpl_cat ~ "Federal Poverty Level",
      highgrade_tvis_cat ~ "Caregiver Education",
      sc_age_years ~ "Child Age",
      sc_race_cat ~ "Child Race",
      sc_sex_cat ~ "Child Sex"
    )) |>
  bold_p(t = 0.05)

# Create a summary table for the Depression model
tbl_dep <- tbl_regression(model_dep, exponentiate = TRUE, label = list(
      ace_counts ~ "ACEs",
      fpl_cat ~ "Federal Poverty Level",
      highgrade_tvis_cat ~ "Caregiver Education",
      sc_age_years ~ "Child Age",
      sc_race_cat ~ "Child Race",
      sc_sex_cat ~ "Child Sex"
    )) |>
  bold_p(t = 0.05)

# Create a summary table for the Behavioral Issues model
tbl_beh <- tbl_regression(model_beh, exponentiate = TRUE, label = list(
      ace_counts ~ "ACEs",
      fpl_cat ~ "Federal Poverty Level",
      highgrade_tvis_cat ~ "Caregiver Education",
      sc_age_years ~ "Child Age",
      sc_race_cat ~ "Child Race",
      sc_sex_cat ~ "Child Sex"
    )) |>
  bold_p(t = 0.05)

# Merge all tables into one
tbl_merged <- tbl_merge(
  tbls = list(tbl_anx, tbl_dep, tbl_beh),
  tab_spanner = c("Anixety Model", "Depression Model", "Behavioral Issues Model")
)

# Display the combined table
tbl_merged
Characteristic
Anixety Model
Depression Model
Behavioral Issues Model
OR1 95% CI1 p-value OR1 95% CI1 p-value OR1 95% CI1 p-value
ACEs








    0 — —
— —
— —
    1 1.30 1.10, 1.53 0.002 1.23 1.01, 1.49 0.042 1.14 1.00, 1.31 0.054
    2-3 1.49 1.26, 1.76 <0.001 1.45 1.20, 1.74 <0.001 1.11 0.97, 1.27 0.13
    4+ 1.86 1.53, 2.26 <0.001 1.78 1.46, 2.18 <0.001 1.37 1.18, 1.59 <0.001
Child Race








    Non-Hispanic White — —
— —
— —
    Non-Hispanic Black or African American 0.77 0.59, 1.02 0.057 0.94 0.72, 1.25 0.7 1.03 0.86, 1.24 0.7
    Non-Hispanic American Indian or Alaska Native 0.52 0.31, 0.93 0.019 1.35 0.68, 3.06 0.4 0.91 0.57, 1.52 0.7
    Non-Hispanic Asian, Native Hawaiian, or Pacific Islander 0.63 0.43, 0.95 0.021 1.34 0.83, 2.30 0.3 1.14 0.84, 1.58 0.4
    Non-Hispanic Other or mixed race 0.89 0.71, 1.14 0.4 1.01 0.80, 1.29 >0.9 0.93 0.78, 1.10 0.4
    Hispanic or Latino of any race 0.71 0.59, 0.85 <0.001 0.84 0.69, 1.02 0.070 0.90 0.78, 1.05 0.2
Child Sex








    Male — —
— —
— —
    Female 1.34 1.18, 1.51 <0.001 1.42 1.25, 1.62 <0.001 1.21 1.09, 1.35 <0.001
Federal Poverty Level








    Less than 100% — —
— —
— —
    100% to 199% 0.94 0.75, 1.17 0.6 0.89 0.71, 1.12 0.3 0.91 0.77, 1.08 0.3
    200% to 299% 1.03 0.82, 1.31 0.8 0.84 0.66, 1.06 0.14 0.83 0.70, 0.99 0.043
    300% to 399% 0.82 0.64, 1.04 0.10 0.79 0.61, 1.01 0.062 0.68 0.57, 0.82 <0.001
    400% or greater 0.91 0.72, 1.14 0.4 0.69 0.54, 0.87 0.002 0.71 0.59, 0.84 <0.001
Caregiver Education








    Less than high school — —
— —
— —
    High school (including vocational) 1.18 0.77, 1.77 0.4 1.10 0.71, 1.64 0.7 1.01 0.73, 1.37 >0.9
    Some college or associate degree 1.18 0.77, 1.74 0.4 1.13 0.74, 1.68 0.6 0.99 0.72, 1.34 >0.9
    College degree or higher 1.30 0.85, 1.93 0.2 1.03 0.67, 1.53 >0.9 0.87 0.63, 1.17 0.4
Child Age 0.97 0.95, 0.99 <0.001 0.96 0.93, 0.98 <0.001 0.86 0.84, 0.87 <0.001
1 OR = Odds Ratio, CI = Confidence Interval
Expand Code 
# Combine the plots into a list vector
all.models <- list()
all.models[[1]] <- model_anx
all.models[[2]] <- model_dep
all.models[[3]] <- model_beh

# Plot the combined version
plot_models(all.models, axis.lim = c(0.2, 5), dot.size = 1, show.p = TRUE, ci.lvl = 0.95, title = "Odds Ratios by Model", legend.title = "Model", m.labels = c("Anxiety", "Depression", "Behavioral Issues"))

Expand Code 
# Logistic regression model with CFI_dich as a moderator for Anxiety
model_anx_mod <- glm(Anxiety ~ ace_counts * CFI_dich + sc_race_cat + sc_sex_cat + fpl_cat + highgrade_tvis_cat + sc_age_years, data = neurodivergent_data, family = binomial)

# Create a summary table for the Anixety model
tbl_anx_mod <- tbl_regression(model_anx_mod, exponentiate = FALSE, label = list(
      ace_counts ~ "ACEs",
      CFI_dich ~ "Child Flourishing Index",
      fpl_cat ~ "Federal Poverty Level",
      highgrade_tvis_cat ~ "Caregiver Education",
      sc_age_years ~ "Child Age",
      sc_race_cat ~ "Child Race",
      sc_sex_cat ~ "Child Sex"
    )) |>
  bold_p(t = 0.05)

# Logistic regression model with CFI_dich as a moderator for Depression
model_dep_mod <- glm(Depression ~ ace_counts * CFI_dich + sc_race_cat + sc_sex_cat + fpl_cat + highgrade_tvis_cat + sc_age_years, data = neurodivergent_data, family = binomial)

# Create a summary table for the Depression model
tbl_dep_mod <- tbl_regression(model_dep_mod, exponentiate = FALSE, label = list(
      ace_counts ~ "ACEs",
      CFI_dich ~ "Child Flourishing Index",
      fpl_cat ~ "Federal Poverty Level",
      highgrade_tvis_cat ~ "Caregiver Education",
      sc_age_years ~ "Child Age",
      sc_race_cat ~ "Child Race",
      sc_sex_cat ~ "Child Sex"
    )) |>
  bold_p(t = 0.05)

# Logistic regression model with CFI_dich as a moderator for Behavioral Issues
model_beh_mod <- glm(Behavioral_Issues ~ ace_counts * CFI_dich + sc_race_cat + sc_sex_cat + fpl_cat + highgrade_tvis_cat + sc_age_years, data = neurodivergent_data, family = binomial)

# Create a summary table for the Behavioral Issues model
tbl_beh_mod <- tbl_regression(model_beh_mod, exponentiate = FALSE, label = list( 
      ace_counts ~ "ACEs",
      CFI_dich ~ "Child Flourishing Index",
      fpl_cat ~ "Federal Poverty Level",
      highgrade_tvis_cat ~ "Caregiver Education",
      sc_age_years ~ "Child Age",
      sc_race_cat ~ "Child Race",
      sc_sex_cat ~ "Child Sex"
    )) |>
  bold_p(t = 0.05)

# Merge all tables into one
mod_tbl_merged <- tbl_merge(
  tbls = list(tbl_anx_mod, tbl_dep_mod, tbl_beh_mod),
  tab_spanner = c("Anixety Model", "Depression Model", "Behavioral Issues Model")
)

# Display the combined table
mod_tbl_merged
1
The interaction anxiety model did not yield significant findings, so we do not proceed with post-hoc analyses.
2
The interaction depression model did yield significant findinds, so we do proceed with post-hoc analyses.
3
The interaction behavioral issues model did yield significant findinds, so we do proceed with post-hoc analyses.
Characteristic
Anixety Model
Depression Model
Behavioral Issues Model
log(OR)1 95% CI1 p-value log(OR)1 95% CI1 p-value log(OR)1 95% CI1 p-value
ACEs








    0 — —
— —
— —
    1 0.26 0.09, 0.43 0.003 0.19 -0.01, 0.39 0.065 0.14 0.00, 0.28 0.043
    2-3 0.38 0.20, 0.55 <0.001 0.32 0.13, 0.51 <0.001 0.11 -0.03, 0.25 0.12
    4+ 0.58 0.38, 0.78 <0.001 0.54 0.33, 0.74 <0.001 0.31 0.16, 0.46 <0.001
Child Flourishing Index








    Not Flourishing — —
— —
— —
    Flourishing -0.99 -1.4, -0.50 <0.001 -1.4 -2.2, -0.61 <0.001 -0.71 -1.3, -0.13 0.013
Child Race








    Non-Hispanic White — —
— —
— —
    Non-Hispanic Black or African American -0.27 -0.54, 0.01 0.050 -0.06 -0.33, 0.23 0.7 0.03 -0.15, 0.22 0.7
    Non-Hispanic American Indian or Alaska Native -0.66 -1.2, -0.07 0.018 0.28 -0.40, 1.1 0.5 -0.11 -0.57, 0.41 0.7
    Non-Hispanic Asian, Native Hawaiian, or Pacific Islander -0.48 -0.85, -0.07 0.016 0.29 -0.19, 0.83 0.3 0.11 -0.20, 0.44 0.5
    Non-Hispanic Other or mixed race -0.13 -0.36, 0.12 0.3 0.01 -0.23, 0.25 >0.9 -0.07 -0.24, 0.10 0.4
    Hispanic or Latino of any race -0.34 -0.52, -0.16 <0.001 -0.18 -0.37, 0.02 0.075 -0.10 -0.25, 0.05 0.2
Child Sex








    Male — —
— —
— —
    Female 0.30 0.18, 0.42 <0.001 0.36 0.23, 0.49 <0.001 0.19 0.09, 0.30 <0.001
Federal Poverty Level








    Less than 100% — —
— —
— —
    100% to 199% -0.06 -0.28, 0.16 0.6 -0.11 -0.34, 0.12 0.4 -0.11 -0.29, 0.06 0.2
    200% to 299% 0.03 -0.20, 0.27 0.8 -0.18 -0.42, 0.06 0.14 -0.20 -0.38, -0.02 0.026
    300% to 399% -0.19 -0.43, 0.05 0.12 -0.24 -0.50, 0.01 0.062 -0.40 -0.59, -0.21 <0.001
    400% or greater -0.09 -0.32, 0.14 0.5 -0.37 -0.61, -0.14 0.002 -0.37 -0.55, -0.19 <0.001
Caregiver Education








    Less than high school — —
— —
— —
    High school (including vocational) 0.14 -0.29, 0.55 0.5 0.09 -0.34, 0.50 0.7 -0.01 -0.34, 0.30 >0.9
    Some college or associate degree 0.13 -0.30, 0.52 0.5 0.12 -0.30, 0.51 0.6 -0.04 -0.36, 0.27 0.8
    College degree or higher 0.22 -0.21, 0.61 0.3 0.02 -0.41, 0.42 >0.9 -0.17 -0.49, 0.14 0.3
Child Age -0.03 -0.05, -0.01 0.002 -0.04 -0.07, -0.02 <0.001 -0.15 -0.17, -0.14 <0.001
ACEs * Child Flourishing Index








    1 * Flourishing -0.23 -0.91, 0.46 0.5 -0.13 -1.3, 1.0 0.8 -0.52 -1.3, 0.24 0.2
    2-3 * Flourishing 0.10 -0.58, 0.80 0.8 1.3 0.27, 2.5 0.017 -0.54 -1.3, 0.20 0.2
    4+ * Flourishing 0.54 -0.44, 1.7 0.3 1.0 -0.08, 2.2 0.077 -1.1 -2.1, -0.19 0.019
1 OR = Odds Ratio, CI = Confidence Interval
Expand Code 
# Create a summary table for the Depression model
tbl_dep <- model_dep_mod |> 
  tbl_regression(
    tidy_fun = tidy_margins,
    estimate_fun = scales::label_percent(style_positive = "plus"),
    label = list(
      ace_counts ~ "ACEs",
      CFI_dich ~ "Child Flourishing Index",
      fpl_cat ~ "Federal Poverty Level",
      highgrade_tvis_cat ~ "Caregiver Education",
      sc_age_years ~ "Child Age",
      sc_race_cat ~ "Child Race",
      sc_sex_cat ~ "Child Sex"
    )
  ) |>
  bold_p(t = 0.05) |>  # Bold significant p-values
  modify_header(label = "**Variable**") |>  # Adjust the header for the variable column
  modify_spanning_header(c(estimate, ci, p.value) ~ "**Depression Model**") |>  # Custom header
  bold_labels()  # Bold variable labels

# Create a summary table for the Behavioral Issues model
tbl_beh <- model_beh_mod |> 
  tbl_regression(
    tidy_fun = tidy_margins,
    estimate_fun = scales::label_percent(style_positive = "plus"),
    label = list(
      ace_counts ~ "ACEs",
      CFI_dich ~ "Child Flourishing Index",
      fpl_cat ~ "Federal Poverty Level",
      highgrade_tvis_cat ~ "Caregiver Education",
      sc_age_years ~ "Child Age",
      sc_race_cat ~ "Child Race",
      sc_sex_cat ~ "Child Sex"
    )
  ) |>
  bold_p(t = 0.05) |>  # Bold significant p-values
  modify_header(label = "**Variable**") |>  # Adjust the header for the variable column
  modify_spanning_header(c(estimate, ci, p.value) ~ "**Behavioral Issues Model**") |>  # Custom header
  bold_labels()  # Bold variable labels

# Merge all tables into one
tbl_combined <- tbl_merge(
  tbls = list(tbl_dep, tbl_beh),
  tab_spanner = c("Depression Model", "Behavioral Issues Model")
)

# Display the combined table
tbl_combined
Variable
Depression Model
Behavioral Issues Model
Average Marginal Effects 95% CI1 p-value Average Marginal Effects 95% CI1 p-value
ACEs





    0 — —
— —
    1 +2.871% -0.232%, +5.974% 0.070 +1.579% -0.1018%, +3.2595% 0.066
    2-3 +5.292% +2.399%, +8.184% <0.001 +1.179% -0.4948%, +2.8521% 0.2
    4+ +7.839% +4.881%, +10.797% <0.001 +3.232% +1.4625%, +5.0025% <0.001
Child Flourishing Index





    Not Flourishing — —
— —
    Flourishing -12.262% -19.286%, -5.237% <0.001 -21.544% -27.7703%, -15.3176% <0.001
Federal Poverty Level





    Less than 100% — —
— —
    100% to 199% -1.300% -4.019%, +1.418% 0.3 -1.161% -2.9205%, +0.5992% 0.2
    200% to 299% -2.211% -5.134%, +0.712% 0.14 -2.154% -4.0256%, -0.2828% 0.024
    300% to 399% -3.061% -6.255%, +0.133% 0.060 -4.562% -6.6859%, -2.4373% <0.001
    400% or greater -4.935% -7.915%, -1.954% 0.001 -4.164% -6.0833%, -2.2442% <0.001
Caregiver Education





    Less than high school — —
— —
    College degree or higher +0.246% -5.459%, +5.951% >0.9 -1.967% -5.4119%, +1.4780% 0.3
    High school (including vocational) +1.254% -4.503%, +7.012% 0.7 -0.148% -3.6319%, +3.3368% >0.9
    Some college or associate degree +1.578% -4.058%, +7.214% 0.6 -0.406% -3.8201%, +3.0072% 0.8
Child Age -0.581% -0.901%, -0.262% <0.001 -1.804% -1.9761%, -1.6315% <0.001
Child Race





    Non-Hispanic White — —
— —
    Hispanic or Latino of any race -2.475% -5.305%, +0.355% 0.086 -1.218% -3.0020%, +0.5661% 0.2
    Non-Hispanic American Indian or Alaska Native +3.419% -4.706%, +11.543% 0.4 -1.273% -7.2925%, +4.7462% 0.7
    Non-Hispanic Asian, Native Hawaiian, or Pacific Islander +3.499% -2.070%, +9.068% 0.2 +1.218% -2.1890%, +4.6259% 0.5
    Non-Hispanic Black or African American -0.780% -4.571%, +3.010% 0.7 +0.381% -1.7004%, +2.4631% 0.7
    Non-Hispanic Other or mixed race +0.103% -3.066%, +3.273% >0.9 -0.848% -2.9258%, +1.2294% 0.4
Child Sex





    Male — —
— —
    Female +4.743% +3.087%, +6.399% <0.001 +2.165% +0.9965%, +3.3334% <0.001
1 CI = Confidence Interval
Expand Code 
# Plot marginal effects for the Depression model
plot_model(model_dep_mod, type = "int", title = "Interaction: ace_counts * CFI_dich (Depression)")

Expand Code 
# Plot marginal effects for the Behavioral Issues model
plot_model(model_beh_mod, type = "int", title = "Interaction: ace_counts * CFI_dich (Behavioral Issues)")

Session Information

To enhance reproducibility, details about the working environment used for these analyses can be found below.

Expand Code 
sessionInfo()
R version 4.4.1 (2024-06-14)
Platform: aarch64-apple-darwin20
Running under: macOS Sonoma 14.6.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/Chicago
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] sjmisc_2.8.10        sjlabelled_1.2.0     patchwork_1.2.0     
 [4] broom.helpers_1.17.0 ggstats_0.7.0        gtsummary_2.0.3     
 [7] margins_0.3.28       sjPlot_2.8.16        webshot2_0.1.1      
[10] kableExtra_1.4.0     broom_1.0.6          Hmisc_5.1-2         
[13] knitr_1.48           lubridate_1.9.3      forcats_1.0.0       
[16] stringr_1.5.1        purrr_1.0.2          readr_2.1.5         
[19] tibble_3.2.1         ggplot2_3.5.1        tidyverse_2.0.0     
[22] tidyr_1.3.1          dplyr_1.1.4          pacman_0.5.1        

loaded via a namespace (and not attached):
 [1] gridExtra_2.3       rlang_1.1.4         magrittr_2.0.3     
 [4] snakecase_0.11.1    prediction_0.3.18   compiler_4.4.1     
 [7] systemfonts_1.1.0   vctrs_0.6.5         pkgconfig_2.0.3    
[10] crayon_1.5.3        fastmap_1.2.0       backports_1.5.0    
[13] labeling_0.4.3      effectsize_0.8.9    utf8_1.2.4         
[16] PupillometryR_0.0.5 promises_1.3.0      rmarkdown_2.27     
[19] markdown_1.13       tzdb_0.4.0          haven_2.5.4        
[22] ps_1.7.6            bit_4.0.5           xfun_0.45          
[25] labelled_2.13.0     jsonlite_1.8.8      highr_0.11         
[28] later_1.3.2         ggeffects_1.7.0     parallel_4.4.1     
[31] cluster_2.1.6       R6_2.5.1            RColorBrewer_1.1-3 
[34] stringi_1.8.4       rpart_4.1.23        Rcpp_1.0.12        
[37] parameters_0.22.0   base64enc_0.1-3     nnet_7.3-19        
[40] timechange_0.3.0    tidyselect_1.2.1    rstudioapi_0.16.0  
[43] yaml_2.3.9          websocket_1.4.1     processx_3.8.4     
[46] bayestestR_0.13.2   withr_3.0.0         evaluate_0.24.0    
[49] foreign_0.8-86      xml2_1.3.6          pillar_1.9.0       
[52] checkmate_2.3.1     insight_0.20.1      generics_0.1.3     
[55] vroom_1.6.5         chromote_0.2.0      hms_1.1.3          
[58] commonmark_1.9.1    munsell_0.5.1       scales_1.3.0       
[61] glue_1.8.0          tools_4.4.1         data.table_1.15.4  
[64] grid_4.4.1          cards_0.3.0         datawizard_0.11.0  
[67] colorspace_2.1-0    performance_0.12.0  htmlTable_2.4.2    
[70] Formula_1.2-5       cli_3.6.3           fansi_1.0.6        
[73] viridisLite_0.4.2   svglite_2.1.3       gt_0.11.1          
[76] sjstats_0.19.0      gtable_0.3.5        sass_0.4.9         
[79] digest_0.6.36       htmlwidgets_1.6.4   farver_2.1.2       
[82] htmltools_0.5.8.1   lifecycle_1.0.4     bit64_4.0.5        
[85] MASS_7.3-60.2

This document was made by Adrián Medina as part of the Slopen Laboratory, based on the Lumo quarto template.